initialPh=7.4;

Hmolar[initialPh_]:=(10^(-initialPh));

initialCO2mmHg=40;

co2SolubilityMilliMolarPerMmhg=0.03;

CO2molarFromMillimetresHg[CO2mmHg_]:=((CO2mmHg co2SolubilityMilliMolarPerMmhg)/1000);

CO2mmHgFromMolar[CO2molarParam_]:=1000 (CO2molarParam/co2SolubilityMilliMolarPerMmhg)

pKfromTextbook=6.1;

KfromTextbook=(10^(-pKfromTextbook));

initialCO2molar=CO2molarFromMillimetresHg[initialCO2mmHg];

HCO3molar[CO2molar_,pH_]:=CO2molar (KfromTextbook/(10^(-pH)));

initialHCO3molar=HCO3molar[initialCO2molar,initialPh];

pCO2mmHg[pH_,HCO3Molar_,K_]:=((10^(-pH))/K)((1000HCO3Molar)/co2SolubilityMilliMolarPerMmhg);

initialHmolar=Hmolar[initialPh];

KbyInitial=((initialHmolar initialHCO3molar)/initialCO2molar);

pKbyInitial = -Log10[KbyInitial];

contourPlot1 =ContourPlot[pCO2mmHg[pH,HCO3/1000,KfromTextbook],{pH,6.8,8},{HCO3,0,100},
Contours->{20,40,80}
]

plot3d1=Plot3D[pCO2mmHg[pH,HCO3milliMolar/1000,KfromTextbook],{pH,6.8,8},{HCO3milliMolar,0,100},
Mesh->None,
ColorFunction->Hue,
ImageSize->Full,
PlotStyle->Directive[Opacity[0.8]],
AxesLabel->{pH,HCO3,CO2}
];

J[CO2_,H_,HCO3_,K_]:=
With[
{kb=1,kf=K},
kf CO2-kb H HCO3];

s2[initialHmolar_,initialCO2molar_,initialHCO3molar_,co2source_,co2sink_,protonSource_,timeSpan_]:=NDSolve[
{
H'[t]==J[CO2[t],H[t],HCO3[t],KfromTextbook]+protonSource,
HCO3'[t]==J[CO2[t],H[t],HCO3[t],KfromTextbook],
CO2'[t]==-J[CO2[t],H[t],HCO3[t],KfromTextbook]+co2source-co2sink,
H[0]==initialHmolar,
CO2[0]==initialCO2molar,
HCO3[0]==initialHCO3molar
},
{H,HCO3,CO2},{t,0,timeSpan}];

Manipulate[
With[{s3=s2[Hmolar[initialPh],CO2molarFromMillimetresHg[initialCO2mmHg],initialHCO3mM/1000,co2source,co2sink,protonSource,timeSpan]},
{
Plot[Evaluate[-Log10[H[t]]/.s3],{t,0,timeSpan},PlotRange->{{0,timeSpan},{6.8,8.0}}],
Plot[Evaluate[ 1000CO2[t]/co2SolubilityMilliMolarPerMmhg/.s3],{t,0,timeSpan},PlotRange->{{0,timeSpan},{0,120}}],
Plot[Evaluate[1000 HCO3[t]/.s3],{t,0,timeSpan},PlotRange->{{0,timeSpan},{0,100}}],
ParametricPlot[Evaluate[{-Log10[H[t]],1000 HCO3[t]}/.s3],{t,0,timeSpan},PlotRange->{{6.8,8.0},{0,100}},AspectRatio->1,
ColorFunction->(ColorData["Rainbow"][#3]&),MaxRecursion->0,PlotStyle->Thick],
{
-Log10[H[timeSpan]],
1000HCO3[timeSpan]
}/.s3
}],
{{initialPh,7.4}, 7.0, 7.8},
{{initialCO2mmHg,40}, 20, 80},
{{initialHCO3mM,1000initialHCO3molar}, 0, 100},
{co2source,0,10 (10^(-6))},
{co2sink,0,10 (10^(-6))},
{protonSource,0,2 (10^(-9))},
{timeSpan, 100,1000},
ContentSize->400,
SaveDefinitions->True
]

Manipulate[
With[{s3=s2[Hmolar[p[[1]]],CO2molarFromMillimetresHg[initialCO2mmHg],p[[2]]/1000,co2source,co2sink,protonSource,timeSpan]},

TableForm[{
Show[contourPlot1,ParametricPlot[Evaluate[{-Log10[H[t]],1000 HCO3[t]}/.s3],{t,0,timeSpan},PlotRange->{{6.8,8.0},{0,100}},AspectRatio->1,
ColorFunction->(ColorData["Rainbow"][#3]&),MaxRecursion->0,PlotStyle->Thick]
],
Plot[Evaluate[{-Log10[H[t]]}/.s3],{t,0,timeSpan},PlotRange->{{0,timeSpan},{6.8,8.0}}],
Plot[Evaluate[{  1000CO2[t]/co2SolubilityMilliMolarPerMmhg}/.s3],{t,0,timeSpan},PlotRange->{{0,timeSpan},{0,120}}],
Plot[Evaluate[{1000 HCO3[t]}/.s3],{t,0,timeSpan},PlotRange->{{0,timeSpan},{0,100}}],
{
-Log10[H[timeSpan]],
1000HCO3[timeSpan]
}/.s3,
p
}]
],
{{p, {7.4,23}}, Locator},
{{initialCO2mmHg,40}, 20, 80},
{co2source,0,10 (10^(-6))},
{co2sink,0,10 (10^(-6))},
{protonSource,0,2 (10^(-9))},
{timeSpan, 100,1000},
SaveDefinitions->True
]

parametricPlot[sol_, timeSpan_]:=
ParametricPlot3D[
Evaluate[{-Log10[H[t]],1000HCO3[t],1000CO2[t]/co2SolubilityMilliMolarPerMmhg}/.sol],
{t,0,timeSpan},
ColorFunction->(ColorData["Rainbow"][1-#3]&),
MaxRecursion->0,
PlotStyle->Thick,
BoxRatios->{1,1,1},
PlotRange->{{6.8,8},{0,100},{0,400}},
ImageSize->Full
];

Manipulate[
With[{s3=s2[Hmolar[initialPh],CO2molarFromMillimetresHg[initialCO2mmHg],initialHCO3milliMolar/1000,co2source,co2sink,protonSource,timeSpan]},
Show[plot3d1,parametricPlot[s3,timeSpan]]
],
{{initialPh,7.7}, 6.8, 7.8},
{{initialCO2mmHg,40}, 20, 400},
{{initialHCO3milliMolar,1000initialHCO3molar}, 0, 100},
{co2source,0,10 (10^(-6))},
{co2sink,0,10 (10^(-6))},
{protonSource,0,2 (10^(-9))},
{timeSpan, 100,1000},
ContentSize->400,
SaveDefinitions->True
]